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REMEDIATING NON-POSITIVE DEFINITE STATE COVARIANCES 
FOR COLLISION PROBABILITY ESTIMATION 


Doyle T. Hall,” Matthew D. Hejduk,t and Lauren C. Johnson? 


The NASA Conjunction Assessment Risk Analysis team estimates the probability 
of collision (P.) for a set of Earth-orbiting satellites. The P. estimation software 
processes satellite position+velocity states and their associated covariance matri- 
ces. On occasion, the software encounters non-positive definite (NPD) state co- 
variances, which can adversely affect or prevent the P, estimation process. Inter- 
polation inaccuracies appear to account for the majority of such covariances, alt- 
hough other mechanisms contribute also. This paper investigates the origin of 
NPD state covariance matrices, three different methods for remediating these co- 
variances when and if necessary, and the associated effects on the P, estimation 
process. 


INTRODUCTION 


The NASA Conjunction Assessment Risk Analysis (CARA) team estimates the probability of 
collision (P.) for a specific set of high-value Earth-orbiting satellites. For each conjunction, the 
CARA software system calculates P, using inertial-frame position+velocity states and associated 
covariance matrices for both the primary and secondary satellites, specified at the time of closest 
approach (TCA). Orbit determination (OD) processes! typically provide positive definite (PD) 
covariances, or at least positive semi-definite (PSD) covariance matrices — characterized by real 
eigenvalues that are all > 0 or = 0, respectively. On occasion, however, the CARA system encoun- 
ters non-positive definite (NPD) state covariances — most often characterized by one slightly neg- 
ative eigenvalue. 


NPD covariances can conceivably arise in a variety of ways. For instance, interpolating from 
an ephemeris of covariances to obtain TCA estimates can introduce inaccuracies that yield NPD 
matrices.” Inaccuracies can also be introduced when converting from one state representation sys- 
tem to another (such as between orbital elements and inertial state vectors) using first-order or 
finite-difference transformation approximations. Computational precision limitations can similarly 
create NPD covariances, even when performing otherwise exact transformations (such as basis- 
vector rotations). NPD matrices can also arise from truncation errors that occur when numerically 
storing or transmitting covariances using data structures with a limited bit length or number of 
significant figures. Preliminary analysis indicates that interpolation inaccuracies account for the 
majority of NPD covariances currently encountered by the CARA system, although other mecha- 
nisms may contribute as well. 
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From a purely computational point of view, P. estimation processes do not always require fully 
PD state covariance matrices for both of the objects involved in a conjunction. For instance, when 
the primary and secondary state covariances are uncorrelated, P, values can be estimated using only 
their joint or summed covariances.*?_ The CARA system occasionally processes conjunctions that 
have an NPD primary or secondary state covariance, but when summed yield a PD joint covariance 
matrix — a situation that raises concern but does not actually prevent the P. computation. Similarly, 
the often-used “2D P.” approximation>* only requires PD status for a marginalized, 2x2 relative 
position covariance at TCA, regardless of the status of the original state covariances. Finally, 
Monte Carlo (MC) methods? require repetitively sampling state probability distribution functions 
(PDFs) for both the primary and secondary objects, but these sampling computations only require 
PSD covariances. 


Despite these situations, the CARA system does sometimes encounter conjunctions in which 
NPD covariances prevent the computation of P. estimates or related quantities. For these situations, 
CARA’s current approach is to “repair” the matrix as required for the specific calculation at hand. 
CARA has investigated three remediation methods for this task: Spectrum Shifting entails adding 
a small offset to the entire set of covariance eigenvalues, and then reconstructing a remediated 
covariance matrix using the original eigenvectors. Higham Remediation® employs a method devel- 
oped for the financial industry to remediate NPD correlation matrices, by finding the closest PSD 
matrix in terms of Frobenius norm. This approach requires that the covariance first be transformed 
into a correlation matrix, which has the advantages of being dimensionless and often better condi- 
tioned than the original covariance. Finally, Eigenvalue Clipping sets a minimum limit for the 
eigenvalues of the covariance or correlation matrix, and simply forces smaller eigenvalues to this 
limit. There are different ways of establishing the eigenvalue limit. For instance, remediating NPD 
position covariances can employ physics-based considerations, specifically by using a clipping 
limit based on the precision of the original OD metric observations, or alternatively, a limit that 
corresponds to a small percentage of the satellite’s hard-body radius. 


This paper investigates the frequency that NPD covariance matrices occur by analyzing a large 
number of conjunctions archived in CARA’s database, as well as describing three different methods 
for remediating these covariances when and if necessary, and the associated effects on the P, esti- 
mation process. 


ANALYSIS OF ARCHIVED CONJUNCTIONS 


This section describes the analysis of actual conjunctions contained in CARA’s data archive, to 
determine the frequency and circumstances in which improper NPD covariances occur. 


Conjunction Data Supplied to CARA 


Most conjunctions represent single, isolated close approach events between two satellites. 
These occur when the magnitude of the relative position vector, |r()|, reaches a minimum in time. 
The U.S. Air Force maintains orbital states for a large catalog of trackable satellites, enabling pre- 
dictions of such close approaches for CARA’s set of high-value primary satellites.’ Conjunction 
analysis occurs when any cataloged secondary satellite is predicted to make an incursion into a 
predefined screening volume centered on a primary.”* For each conjunction, the CARA system 
receives and processes raw data similar to that contained in a standard conjunction data message 
(CDM) — produced in a previous and separate analysis. For purposes of this paper, the input data 
will be referred to hereafter as a CDM. Each CDM contains the time of close approach (fea or 
TCA), and many other quantities related to the orbital states and associated uncertainties of the 
primary and secondary objects for the conjunction. 


ECI Frame Satellite States and Covariances 


For instance, the CDM contains data that can be readily converted to the Earth-centered inertial 
(ECT) reference frame’ mean position+velocity states for both the primary and secondary, Xp(tca) 
and X;(fca), respectively, where x(f) = [r(t) v(t)]. CDM’s also contain state covariance data, origi- 
nally expressed in the Cartesian radial, in-track and cross-track (RIC) reference frame” (also some- 
times called the UVW frame). The RIC covariance matrices can be transformed to the ECI state 
covariances”!° for the primary and secondary, denoted here as C,(fca) and C,(tea), respectively, 
which have matrix dimensions of 6x6. These transformed quantities can then be used to define the 
relative ECI state at close approach, X(fca) = Xs(tca) — Xp(tca), and the associated combined covariance 
C(tea) = Cs(tea) + Cp(tea). The CARA software system employs TCA ECI states and covariances to 
estimate collision probabilities. 


Marginalized Position and Velocity Covariances 


The 6x6 ECI covariance matrices discussed can be decomposed into three 3x3 sub-matrices as 
follows 


A B 
Cal f (1) 


where A denotes the marginalized position-position covariance, E the marginalized velocity-ve- 
locity covariance, and B contains the position-velocity cross covariance terms. This report denotes 
covariance matrices and sub-matrices using all upper case letter symbols. 


Correlation Matrices 


Under most circumstances state covariances can be transformed into correlation matrices, de- 
noted here using lower case symbols. These have unit diagonals, and can be calculated as follows: 
the correlation matrix, c, associated with covariance matrix, C, is given by'! 


c=DCD (2) 


where D denotes a diagonal matrix with elements Dj; = ;; (Cii)""”, and 6;; denotes the Kronecker 
delta function. This conversion yields a proper, unit-diagonal correlation matrix only if all of the 
diagonal elements, C;,;, are positive, which may not be true for an approximated covariance. In 
conjunction analysis, correlation matrices have the advantages of being dimensionless, and also 
usually much better conditioned than the original covariance. For instance, a 6x6 ECI state covar- 
lance, as given in equation (1), has elements with dimensions of m? in the A submatrix, ms” in 
the E submatrix, and m’s"! in the B submatrix. The six eigenvalues of such a mixed-dimension 
covariance are difficult to interpret physically, and can also span a much larger numerical range 
than those for the associated correlation matrix. Matrix condition numbers (calculated using the 
cond function of the Matlab software system) for the ECI state covariances C,(tca) and Cs(tca) were 
found to have typical values of order ~3x10! in the archive analysis, but condition numbers for the 
associated correlations matrices ¢,(tca) and €s(tca) were ~3x10’, six orders of magnitude smaller. 
The marginalized position covariances A,(tca) and As(tca), had significantly better typical condition 
numbers of ~3x10°, and their correlation matrices a,(fca) and a(tca) had even better condition num- 
bers ~10°. 


Equinoctial Element States and Covariances 


The ECI states and covariances for the primary and secondary objects can also be converted to 
equinoctial element mean states and covariances” !° denoted here using the primed symbols x'i(tca) 
and C'x(tca), respectively, where k represents either p or s. For this study, the Jacobian matrices 
required for ECIl<>equinoctial covariance transformations were calculated numerically using a fi- 
nite differencing algorithm. The equinoctial state covariances were found to be somewhat better 
conditioned than their associated ECI state covariances, with condition numbers typically factors 
of about 10 to 20 smaller. Equinoctial state correlation matrices were found to be the best condi- 
tioned among all of the 6x6 matrices analyzed here, with typical condition numbers ~3x 107. 


NPD Number and Status Indicators 


Properly formulated and constructed covariance and correlation matrices cannot be NPD.”°"! 
Trying to use NPD matrices without any remediation to calculate P. values could prevent the cal- 
culation or potentially lead to inaccurate results. To prevent this, such NPD matrices must first be 
detected, and this section defines two indicators designed for this task. The first indicator of the 
NPD status of an NxN covariance matrix, C, is given by the NPD number, n(C), an integer quantity 
defined as the summed number of nonpositive eigenvalues. (All eigenvalues in this analysis were 
calculated using Matlab’s eig function.) PD covariances have n(C) = 0, and NPD covariances have 
0 < n(C) < N. The second indicator is the NPD status, b(C) = min[1, n(C)], a binary quantity 
defined so that PD matrices have b(C) = 0, and NPD matrices have b(C) = 1. 


The NPD number indicator for a correlation matrix, c, can be defined similarly, except in the 
rare cases in which the original NxN covariance matrix, C, has one or more diagonal elements less 
than or equal to zero. As explained above, this prevents the calculation of a proper, unit-diagonal 
correlation matrix. In these cases, the NPD number indicator is defined as n(c) = N+L, where L is 
the number of nonpositive diagonal elements of C. Defined this way, PD correlation matrices have 
n(c) = 0, NPD correlations have 0 < n(c) < N, and, in cases where proper correlations cannot be 
calculated, N < n(c) <2N. The binary NPD status indicator for correlation matrices remains defined 
as b(¢c) = min[1, n(c)]. 


NPD Fractions 


The NPD analysis for a large, archived set of conjunctions can be conveniently summarized in 
the form of NPD fractions, which indicate the fraction of conjunctions involving an NPD matrix. 
For instance, for a set of M conjunctions, the fraction that have NPD primary object ECI state 6x6 
covariances can be expressed as Fp = {Xml b(Cmp)]} 7M, where &,, represents a summation over the 
set of M conjunctions, and C,,,, is the primary’s ECI state covariance at TCA for the m" archived 
conjunction. Likewise the NPD fraction for ECI relative state correlations is f= {2n[b(em)] }/M. 
The NPD fractions {F), fy, Fs. fs, Ff, Fb. fp FS, f%} can all be defined similarly (again with upper- 
case symbols corresponding to covariances, lower-case to correlations, unprimed quantities to ECI 
states, primed quantities to equinoctial states, {p,s} subscripts to primary and secondary object 
states, and no subscripts to relative states). 


NPD Analysis of Two Years of Archived Conjunction Data 


This study reports the analysis of 839,383 conjunctions that occurred during a 2-year period 
between 2015-04-01 and 2017-04-06. Of these, 8,874 were eliminated due to primary or secondary 
ECI state covariance matrices that contained all zeros, or that defined implausibly large 1-sigma 
position uncertainties greatly exceeding one earth radius. NPD number and status indicators, as 
defined above, were calculated for all of the remaining M = 830,509 conjunctions. Specifically, 
these indicators were calculated for the 6x6 ECI and equinoctial state covariance and correlation 


matrices, as well as the 3x3 ECI position covariance and correlation matrices, and then used to 
calculate NPD fractions. 
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Figure 1. NPD fractions plotted as a function of time during the 2-year study period. 


Figure 1 shows the NPD fractions for 6x6 state matrices {Fp, fo, Fs, ft, F'>, fp, F's, f%} plotted as 
a function of time during the 2-year study period using 1-week time bins. In Figure 1, the dashed 
blue lines correspond to covariance matrices and the dotted red lines to correlation matrices; the 
left-side panels to ECI states and the right panels to equinoctial states. Specifically, the top left 
panel of Figure 1 shows NPD fractions for the primary ECI state matrices (F’, dashed blue, f, dotted 
red); the bottom left for secondary ECI matrices (F; dashed blue, f; dotted red); top right for primary 
equinoctial state matrices (F'4, dashed blue, f/, dotted red); and bottom right for secondary equinoc- 
tial state matrices (F’4 dashed blue, f% dotted red). The corresponding NPD fractions for the 3x3 
position matrices are all significantly smaller, and are not plotted here; in fact, the fractions { Fp, fy, 
>» fp} were all found to be zero, except for one primary which had a small number of NPD oc- 
currences (as discussed further below). 


Notably, Figure 1 indicates that NPD fractions observed during the first 13 months of the 2-year 
study period were significantly larger than during the last 11 months, with a dramatic and sudden 
decrease occurring just before 2016-05-01. This improvement occurred because, near this time, 
the number of significant figures used to express RIC covariance elements within the CDMs ex- 
panded from 4 up to 16 — an increased precision that markedly reduced the number of NPD 6x6 
matrices. This demonstrates explicitly that NPD matrices can arise when numerically storing or 
transmitting covariances using data structures with a limited precision due to a limited number of 
significant figures or bit length. 


Because the earlier, lower precision data does not provide a good representation for future con- 
junctions, the remainder of the archive analysis concentrates on the later 11-month period. 


NPD Analysis of High-Precision Archived Conjunction Data 


A total of 433,768 conjunctions occurred during the 11-month period between 2016-05-01 and 
2017-04-06. Of these, 5,179 were eliminated, as before, due to all-zero covariances or implausibly 
large position uncertainties. The remaining M = 428,589 conjunctions comprised 81 unique pri- 
mary satellites, and occurred at a rate of about 8,000 per week. Again, NPD number and status 
indicators were calculated for all of these conjunctions. Table 1 along with Figures 2 and 3 sum- 
marize the results in the form of the NPD fractions {F), fp, Fs, fs, Fb, f'» F %,fs}, as described earlier. 


Table 1. NPD fractions for conjunction covariance and correlation matrices (%). 


Number Num. NPD Fractions (%) for 6x6 State Covariances (F) and Correlations (f) 
F Of : of Cartesian ECI States Equinoctial States 
Primaries Events 
M Fp lp F; Ss F f Fy fo F's | fs 
All 81 428589 6.81 6.76 1.62 0.20 0.12 0.06 6.76 6.76 0.20 0.20 


72 with en <0.01 308971 | 8.4e-2 4.1e-2 1.36 0.14 84e-2 2.0e-2 4.1le-2 41e-2 0.14 0.14 

9withe,>0.68 119618 24.19 24.12 2.28 0.36 0.20 0.17 | 24.12 | 24.12 0.36 0.36 
NPD Fractions (%) for 3x3 Position Covariances (F) and Correlations (f) 

All 81 428589 6.8e-3 6.8e-3 1.2e-3 1.2e-3 1.9e-3 1.9e-3 NA* NA NA. NA 

72 with e,<10? 308971 0 0 6.5e-4 6.5e4 65e-4 65e4 NA NA NA. NA 

9 with ep>0.68 119618 | 2.4e-2 | 2.4e-2  2.5e-3 2.5e-3  5.0e-3  5.0e-3 NA NA NA. NA 


*NA = Not applicable; 3x3 marginalized matrices are only analyzed for Cartesian states. 


The top part of Table 1 lists the results for analyses of 6x6 state matrices, and the bottom part 
summarizes the results for the marginalized 3x3 ECI position matrices. Among the 81 primary 
satellites, 72 were found to have low eccentricity orbits, with e, < 10; the remaining 9 had very 
eccentric orbits, with 0.68 < e, < 0.84. A total of 29,201 of the 428,589 ECI primary state 6x6 
covariances were found to be NPD, representing a fraction of F;, = 6.81%, as listed in the top data 
row of Table 1. The fraction of NPD primary ECI correlation matrices was slightly smaller at f, = 
6.76%, listed adjacently. The bottom part of Table 1 lists the NPD fractions for the marginalized 
3x3 ECI position covariances and correlations {Amp, Amp, Ams, &np, Am, An}, Which are signifi- 
cantly smaller than the NPD fractions for the 6x6 matrices (as discussed further below). 


Figures 2 and 3 show NPD fractions for the 6x6 state matrices plotted as a function of time 
during the 11-month study period. Figure 2 shows the fractions for both primary and secondary 
objects for all of the conjunctions combined, in the same format used in Figure 1 above. Figure 3 
separates the fractions for the primaries into low- and high-eccentricity groups. Specifically, the 
top panels of Figure 3 show NPD fractions for the 72 primaries with e, < 107, and the bottom panels 
for the 9 other primaries with 0.68 < e, < 0.84. Again, the corresponding NPD fractions for the 3x3 
position matrices are all zero, except for one primary (as discussed further in the next section), and 
are not plotted here. 
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Figure 3. 


Empirical Observations from the NPD Analysis 


The following observations can be made from the 11-month archive analysis summarized in 


Table 1 as well as Figures 2 and 3: 


1. 


The 3x3 ECI position covariance and correlation matrices {Aim,p, Amp, Ams, &mp, Am, @n} were 
found to have significantly fewer NPD occurrences than their 6x6 counterparts. In fact, only 
one primary — Themis E (SCN 30798), in a highly eccentric orbit with e, = 0.83 — had any 
NPD 3x3 matrices: specifically, in 29 out of its 8,905 conjunctions. All other primaries had 
3x3 fractions F, =f, = 0 exactly. When including all primaries, the NPD fractions {F,, f,} for 
the 3x3 matrices were both less than 0.01%, but for the 6x6 matrices both approached 7%. In 
other words, the NPD frequency for 3x3 matrices could be considered to be negligibly small 
compared to that of 6x6 matrices. For this reason, the remainder of this section restricts dis- 
cussion to the 6x6 matrices, unless otherwise stated. 

NPD fractions for high-eccentricity objects were found to significantly exceed those for low- 
eccentricity objects (see Figure 3). For example, during the 11-month study period, 23 of the 
72 primaries with e, < 10° were found to have zero NPD fractions for their ECI covariances, 
and only three had F,, > 2%. However, the 9 high-eccentricity primaries all had much larger F, 
values. Analysis shows that inaccurate interpolation of covariance ephemerides generate a 
significant portion of these NPD covariances. Specifically, the ephemerides can have relatively 
large angular spacing near orbital perigee, which degrades the accuracy of the covariance in- 
terpolation. 

Preliminary analysis indicates that primary satellites engaged in maneuvers may also have el- 
evated NPD frequencies, possibly because such objects often require special, customized orbit 
determination processing. Because many CARA primary satellites regularly maneuver, this 
may explain some of the elevated NPD fractions observed for low eccentricity primaries. 

The overall NPD fractions for primary objects typically exceeded those for secondaries by a 
factor of ~2 or more. This is likely related to the fact that CARA primaries, on average, have 
higher eccentricities and maneuver more often than the secondaries. 

Converting ECI state covariances into correlation matrices maintained or reduced the NPD 
fractions for all objects, 1.e., fi < Fi for all k € {p,s}. The ECI covariances {Cnp, Cs, Amp, 
Aims} were found to have all positive diagonal elements without exception for the entire data 
set, allowing calculation of their associated correlation matrices. These ECI covariance—cor- 
relation conversions did not ever produce an NPD correlation from a PD covariance, meaning 
that b(¢mx) < b(Cm) and b(amx) < b(Amx«) for all m € {1...M} andk ¢€ {p,s}. 

Converting ECI state covariances into equinoctial covariances always maintained or reduced 
the NPD fractions, i.e., F; < F% for all k € {p,s}. However, 116 (0.03%) and 79 (0.02%) of the 
equinoctial covariances C’,,, and C’,,; were found to have at least one nonpositive diagonal 
element, respectively, possibly caused by inaccuracies in the ECI—equinoctial transfor- 
mations, which employed Jacobians approximated using a finite differencing scheme. 

The NPD status indicators for the entire set of ECI state correlation matrices and equinoctial 
covariance matrices were found to be identical, and equal to or less than those for the ECI 
covariances, 1.e., b(€m) = b(C'mx) = b(e'mn) < b(Cmx) for all m € {1...M} andk € {p,s}. This 
indicates that the conversions Cin,» > Cmp and Cin,p > C'm,p each had the net effect of maintain- 
ing or reducing NPD occurrences. Furthermore, these conversions never created NPD matrices 
when there was none originally. These NPD status improvements may be related to the better 
condition numbers of the ¢,,, and C’,,, matrices compared to the C,,,» matrices, as discussed 
previously. 


8. The process of summing the primary and secondary ECI covariances, 1.e., Cn = Cis + Cm,p and 
Am = Ams + Amp, did not ever create an NPD combined matrix when there was none originally. 
Also, the NPD status indicators for the relative state matrices were always less than or equal to 
the maximum NPD status for the primary or secondary matrices, i.e., b(Cn) < max[b(Cnp), 
b(Cns)], b(em) < max[b(enmp), b(Cms)], b(Am) < max[b(Amp), b(Ams)], bam) < max[b(amp), 
b(am,s)] for all m € {1...M}. In fact, the covariance summing process often resulted in a PD 
combined covariance for conjunctions in which the primary or secondary covariances were 
NPD, thereby reducing the overall NPD fractions so that F < min[F;, Fs] and f< min[/, fs]. 

9. Although not shown here, NPD frequencies for covariances expressed in the Cartesian radial, 
in-track and cross-track (RIC) reference frame? (also sometimes called the UVW frame) show 
qualitatively the same patterns as those expressed in the Cartesian ECI reference frame. In 
fact, the covariances currently delivered to CARA are expressed in this frame, and were origi- 
nally produced by interpolating between ephemeris times that bracket the conjunction’s TCA 
(c.f., reference 12). 


METHODS OF NPD COVARIANCE REMEDIATION 


The CARA system does sometimes encounter conjunctions in which NPD covariances prevent 
the computation of P, estimates. For instance, the often-used “2D P.” approximation** requires 
PD status for a marginalized, 2x2 relative position covariance at TCA’. Notably, among all of the 
>800,000 conjunctions analyzed here, only one such NPD 2x2 covariance was found, preventing 
2D P, estimation; this involved a high eccentricity primary with e, = 0.68. CARA’s “3D P.” ap- 
proximation’? experiences more frequent NPD issues, because it employs an ephemeris of states 
and covariances spanning the conjunction, and requires PD relative position covariances, A(f), at 
all ephemeris times. Similarly, the relative position Mahalanobis distance at TCA, Mp(tca), which 
can be used to screen efficiently for negligibly small P. values,'* requires PD relative position co- 
variances, A(tca). Finally, CARA’s Monte Carlo P, estimation algorithm requires equinoctial state 
covariances with eigenvalues = 0. 


When the CARA system encounters such problematic cases, the current approach is to remedi- 
ate the NPD matrix to a sufficient level that allows the calculation at hand to be performed. Three 
remediation methods have been investigated for this task: Spectrum Shifting, Higham Remediation, 
and Eigenvalue Clipping. 


Spectrum Shifting Method 


Spectrum Shifting entails adding an offset to the entire set or spectrum of matrix eigenvalues, 
and then constructing a remediated matrix using the original eigenvectors. The first step of this 
NPD remediation process entails describing an NxN matrix, C, (which is assumed here to be real 
and symmetric) using an eigen-decomposition 


C=VAV' (3) 


where V is a unitary matrix comprising the N orthogonal eigenvectors, and A a diagonal matrix of 
the N associated eigenvalues, Ajj = 5;jA;. In this analysis, eigenvalues are sorted into increasing 
order so that A; < Ai+1, and associated eigenvectors are denoted V;. The matrix C requires no reme- 
diation if all of its eigenvalues are positive. If not, the second step defines a shifted set of eigen- 
values, y;, constructed by adding an offset so all become nonnegative: yi = Ai + (Amin| + €), where 
Amin = min[Ai] and ¢ = 0. The third step entails using these shifted eigenvalues to construct a reme- 
diated matrix 


Coon = V EV" © 


where Ii; = 6,7. For ¢ = 0, this remediated matrix is PSD by construction, because Ymin = 0 iden- 
tically. For ¢ > 0, Crem is PD by construction with minimum eigenvalue Ymin = &. CARA’s imple- 
mentation of the spectrum shifting remediation algorithm employs a default value for ¢ that is small 
compared to typical eigenvalues but significantly larger than the computing precision limit. How- 
ever, € can be made even larger if needed (for instance, to ensure that Matlab’s Cholesky decom- 
position function chol evaluates Cyem as PD, in addition to Matlab’s eig function). 


The principal advantage of the spectrum shifting method is simplicity. One disadvantage is that 
the method shifts the entire spectrum of eigenvalues upward, even though usually only one of these 
estimated quantities is originally negative. Another is that the method assumes that the original 
eigenvectors can be used without modification to construct the remediated matrix in equation (4). 
A final disadvantage is that there is no objective, physics-based means of establishing what the size 
of ¢ should be in cases that require a fully PD remediated covariance. 


Higham’s Nearest Correlation Matrix Method 


Higham Remediation employs methods developed for the financial industry to remediate NPD 
correlation matrices® and covariance matrices!!'* by finding the closest PSD matrix in terms of 
Frobenius norm. The CARA software implementation employs the correlation matrix remediation 
algorithm’ as the primary method. So when remediating an approximated NxN covariance matrix, 
C, the first step of the process entails conversion into the associated correlation matrix, c, using 
equation (2). As mentioned previously, this requires that the diagonal elements of C all must be 
positive to yield a proper, unit-diagonal ¢ matrix. In the rare cases where any C;; < 0, Higham’s 
alternate but closely related method" is used as the backup algorithm directly for covariance reme- 
diation. Given a correlation matrix, c, Higham’s algorithm® employs nonlinear optimization to 
estimate the unique PSD matrix with the closest Frobenius matrix norm, |le||? = 2;,[(ci,)*]. The 
resulting remediated correlation matrix, Cem, is proper in that it has unit-diagonal and nonnegative 
eigenvalues. The covariance associated with this remediated correlation can be derived by invert- 
ing equation (2), Crem =D" Crem D'. The conjunction archive analysis indicates that this method 
yields C,em matrices that also usually have nonnegative eigenvalues. However, when that is not the 
case, then again Higham’s alternate method" can be used as a backup algorithm. 


The principal advantages of the Higham remediation methods are that they have been derived 
with mathematical rigor, and multiple versions of the implemented algorithms are posted on-line. 
The disadvantages are that these algorithms are relatively computationally intensive, and designed 
only to estimate PSD matrices — again leaving unaddressed situations requiring a fully PD reme- 
diated matrix. Furthermore, finding the closest matrix in terms of the closest Frobenius norm, while 
mathematically well defined, is not a physics-based criterion. 


Figenvalue Clipping Method 


Eigenvalue Clipping sets a minimum limit for the eigenvalues of a covariance or correlation 
matrix, and simply forces smaller eigenvalues to this limit. The CARA analysis team developed 
this method in part to address situations requiring fully PD covariances, such as 3D P. or Mahalano- 
bis distance estimation. The first step of this remediation process entails eigen-decomposition, as 
in equation (3). Next, the matrix is deemed to need remediation if any of the eigenvalues are found 
to be smaller than a nonnegative eigenvalue clipping limit, Acip, which can be constrained to have 
a sensible size using physics-based considerations as discussed below. Then, a new set of clipped 
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eigenvalues can be defined as y; = max[Aji, Aciip], and the remediated matrix calculated using equa- 
tion (4). To produce PSD remediated matrices Aciip = 0; for fully PD remediation Aci > 0. When 
remediating correlation matrices, a final step employing a transformation in the form of equation 
(2) can be used to ensure that the final result has proper unit diagonal form. 


The principal advantages of eigenvalue clipping are that it is the simplest to implement among 
the three methods, and the physics of hard-body collisions can be used as a basis for remediating 
position covariances. A disadvantage (shared with spectrum shifting) is that this method assumes 
that the original eigenvectors can be reused without modification in equation (4). 


Eigenvalue Clipping for Marginalized Position Covariances 


Physics-based considerations can provide sensible constraints on the eigenvalue clipping limit, 
Actip, for marginalized position covariances. This is because the eigenvalues for such covariances 
can be interpreted? as the squares of the 1-sigma widths of the position uncertainty PDFs, specifi- 
cally 4;= 67. For instance, if a primary object has a 3x3 position covariance, A,, with all positive 
eigenvalues, then the associated PDF can be visualized as a 1-sigma uncertainty ellipsoid with 
principal axes aligned with the eigenvectors (V,1, Vp2, Vp3), and half-width dimensions of (6,1, 
Op,2, Op,3). However, if A, is PSD with 1,1 = 0 and 0 < Ap,2 < Ap3, then the PDF becomes confined 
to a plane, and the 1-sigma surface an elliptical area on that plane. This PSD covariance unrealis- 
tically indicates that there is no uncertainty at all in the object’s position vector along the direction 
of the first eigenvector. NPD covariances are even more unrealistic. They lack realism because 
satellite positions are estimated quantities, based on measurements limited in both number and ac- 
curacy, so actual sigma values must always be real and finite. These unrealistic situations occur 
because A, and its eigenvalues suffer from inaccuracies, which can arise a variety of ways as pre- 
viously discussed. 


Eigenvalue clipping constrains the remediated position covariance to have a 1-sigma uncertainty 
ellipsoid with 6; > (Aciip)'.. One conceivable method of establishing this clipping limit would be 
simply to impose a limit on how accurate future satellite positions could ever be predicted from the 
orbital determination estimation process. For instance, detailed analysis and/or long-term experi- 
ence might allow one to conclude that, for a given set of tracking sensors and measurements, 1- 
sigma uncertainties on future satellite positions could never be reduced below a lower limit, Gjim. 
In this case one could sensibly use a clipping limit of Actip = (Glim)*. Unfortunately, such ovim values 
will likely vary significantly over time and from satellite to satellite; in addition, they may require 
a prohibitively detailed analysis of the original satellite tracking data to estimate confidently. 


Fortunately, considerations based on the physics of collisions can also be used as the basis for 
covariance remediation in these situations. The P, estimation process idealizes collisions as occur- 
ring when the hard-body radii of the primary and secondary objects begin to overlap as they ap- 
proach one another.*° If any of the primary’s estimated eigenvalues are significantly less than the 
square of its hard-body radius, i.e., 0 < A,i << h,’, then random variations in the position vector 
along the eigenvector direction V,,; are very small compared to the object size. As explained in 
more detail below, this means that the probability of collision becomes insensitive to the exact 
value of A), as long as it is significantly smaller than h,”. In other words, all eigenvalues satisfying 
0 < Api << hy” should produce approximately the same P. value. This concept also applies to sec- 
ondary objects. 


To envision this concept, imagine a conjunction with one small primary position eigenvalue, 
ie., 0 < Api << h,’, but no other similarly small eigenvalues i); or As,;.. Next, consider the family 
of related conjunctions in which 4,1 varies over the domain 0 < Api < A)”, but all other conjunction 
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parameters remain fixed. Collision probabilities for this family vary with 4,1. Let P<. denote the 
probability for the family member with A,,; = 0. In this case, although there are no random position 
variations along V,1, the object itself extends along this direction by +h,. For other members of 
this family with A,,1 << h,’, this situation remains essentially unchanged: the object itself extends 
along V,,1 by th, with very small random variations. So these cases should also have collision 
probabilities approximately equal to P..o. Figure 4 demonstrates this explicitly by plotting 2D P. 
values‘ as a function of the ratio 6,,1/hp, for the following conjunction TCA parameters (all speci- 
fied using length units in m, and time units in s): hp = 10, r, = [0 0 0]’, v, = [1070 0]’, V, = [29 
], Op2 = 20, op3 = 100, hs = 1, rs =[1 0 0]’, vs = [0 10* OJ’, Vs = [¥ X Z], 05,1 = 20, 62 = 50, and 
6,3 = 200, where x, ¥, and Z denote standard unit column vectors. For this notional conjunction, 
P.o = 0.0179, and P. values depart from this by less than 0.2% over the range 0 < 6p,1/hp < 0.1 (see 
Figure 4). Analysis indicates that this same pattern holds for many different conjunctions with 
widely varying encounter geometries, both with fabricated parameters (as in this example), or those 
taken directly from CARA’s archive. 


This concept can be exploited for NPD covariance remediation by extending it to the relatively 
small negative eigenvalues that occasionally occur in conjunction analysis. Again, when the CARA 
system encounters NPD covariances, most commonly only one of the eigenvalues is negative, and 
only very slightly so relative to the largest positive eigenvalue. Clipping these slightly negative 
eigenvalues is appropriate if one assumes they represent estimation inaccuracies (such as interpo- 
lation errors), and that a better estimation procedure would have yielded small positive values. 
CARA’s remediation algorithm for 3x3 position covariances employs a clipping limit that corre- 
sponds to a small fraction of the object’s hard-body radius, i.e., Aciip = (Afciip)”, with the clipping 
fraction set to feip = 10 by default, although values over the range 10° < fczip < 10°? produce essen- 
tially equivalent P, estimates in all cases examined. Relative position covariances, A = A, + As, 
can be remediated using the same approach but using the combined hard-body radius, h = hy + hs. 
Notably, this eigenvalue clipping approach produces PD remediated matrices, enabling computa- 
tion of Mahalanobis distances and 3D P- values’? which both require invertible covariances. It also 
enables 2D P, computation in the extremely rare conjunctions where the 2x2 marginalized relative 
position covariance used within that calculation*“ is found to be NPD. 
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Figure 4. P, plotted as a function of the smallest primary covariance sigma value. 
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Eigenvalue Clipping for Monte Carlo Sampling 


CARA’s Monte Carlo P. estimation algorithm requires repetitively sampling the equinoctial 
state probability distribution functions of the primary and secondary objects. Given an object’s 
mean equinoctial state, x’, and covariance, C’, a sampled state, y’, can be written as 


N 
y'=> ov] (5) 
i=] 


where V’; denotes the eigenvectors, 4’; the eigenvalues, and o; normally-distributed random variates 
(generated using Matlab’s randn function). For NPD covariances, this equation will produce sam- 
pled states with imaginary values. However, real-valued sampled states can be approximated using 
an eigenvalue clipping limit of zero, Aciip = 0, as follows 


y-SbyvmatZa 6 
i=l 


The very slight alteration between equations (5) and (6) indicates how simple the software imple- 
mentation of eigenvalue clipping can be, especially compared to Higham remediation. 


Again, most commonly only one of the NPD covariance eigenvalues is negative, and only 
slightly so. In these cases, using equation (6) restricts the distribution of sampled states to a five 
dimensional locus within the six dimensional space of equinoctial states. Similar concepts also 
apply to sampling Monte Carlo states using the equinoctial state correlation matrix c’ rather than 
the covariance C’. This approach may have advantages because of the dimensionless nature and 
better conditioning of the c’ matrices, and is currently being investigated by the CARA team. 


EFFECTS OF COVARIANCE REMEDIATION ON COLLISION PROBABILITIES 


As mentioned previously, from a purely computational point of view, P. estimation processes 
do not always require fully PD state covariance matrices for both of the objects involved in a con- 
junction. This is especially true for the often-used 2D P. approximation>* which only requires PD 
status for a marginalized 2x2 relative position covariance at TCA, regardless of the status of the 
original state covariances. Of the 830,509 archived events analyzed here, only one conjunction 
produced such an NPD 2x2 marginalized covariance that prevented the 2D P, computation. So the 
vast majority of events that have NPD 6x6 covariances for their TCA primary, secondary and/or 
relative states can still be processed to provide a 2D P, estimate without any covariance remediation 
at all. Even the more demanding 3D P, computation’? can be performed in many of these situations 
without invoking any remediation for similar reasons. 


These situations raise concern because of the possibility that the calculated 2D and 3D P, esti- 
mates may be inaccurate because the original covariances were compromised. One way to inves- 
tigate this possibility is to remediate any NPD 6x6 covariances before they are processed for P- 
estimation, and then compare the resulting “P.,em” values to the “Pe unrem” Values calculated without 
any remediation performed. This study uses the 2D P, estimation method to perform such a com- 
parison for the 11-month archive data set described above, comprising 428,589 conjunctions with 
high-precision covariances. When the subset of conjunctions with NPD 6x6 ECI state covariances 
were remediated using eigenvalue clipping with Acip = 0, differences between Porem and Peunrem 
were all found to be less than 0.2% when restricted to the 2,866 conjunctions with Pemax = 
max[P.,rem, Pc.unrem] Values of 10°’ or more. Spectrum shifting remediation produced comparable or 
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smaller differences. Higham remediation also produced comparable differences, except in one 
conjunction (out of 2,866) in which Porem and Peunrem differed by about 1.8%. This demonstrates 
that remediating the full 6x6 ECI covariances before the 2D P, calculation does not substantially 
change P, estimates, regardless of which of the three remediation methods are used. Similar results 
are also found using the 3D P, calculation. 


Other similar tests could also be performed. For instance, instead of directly remediating any 
NPD 6x6 ECI state covariances before the P,. calculation, as done above, the covariances could 
first be converted to equinoctial covariances, remediated in that form, and then converted back into 
ECI state covariances to be used for the “Pe rem” estimates. Similarly, one could first remediate the 
RIC covariances, followed by RIC—>ECI conversions, and then the P, calculation. Converting to 
correlation matrices as an intermediate step could also be investigated. However, as mentioned 
previously, in most cases such remediation is not required from a purely computational point of 
view. CARA’s current approach entails performing remediation only when required by a specific 
calculation. So to calculate a 2D P, estimate, remediation would be performed only for the 2x2 
marginalized covariance required for that specific calculation (i.e., only once for all of the >800,000 
conjunctions analyzed here). 


CONCLUSIONS 


This analysis indicates that the frequency of NPD 6x6 state covariance matrices in CARA’s 
archive decreased dramatically near 2016-05-01, when the number of CDM significant figures was 
increased to improve numerical precision. After this time, the fraction of conjunctions involving 
NPD ECT state covariances for the primary objects decreased to less than 0.1% for low-eccentricity 
primaries. However, almost 25% of the 6x6 ECI covariances for high-eccentricity primaries were 
found to be NPD, likely due at least in part to covariance interpolation inaccuracies. The frequency 
of NPD 3x3 position matrices was significantly smaller than that of 6x6 state matrices. 


Converting covariance state matrices to correlation matrices maintains or decreases the fre- 
quency of NPD occurrences. Converting Cartesian frame state covariances into equinoctial state 
covariances similarly decreases NPD occurrences. 


Collision probability calculations do not always require fully PD state covariance matrices for 
both of the objects involved in a conjunction, because these computations often employ marginal- 
ized covariances of reduced dimension. Of the >800,000 archived conjunctions analyzed here, 
only one produced an NPD marginalized 2x2 covariance required for the 2D P, computation. NPD 
covariances occasionally need to be remediated when calculating Mahalanobis distances, using the 
3D P. approximation, or when performing Monte Carlo simulations to estimate P. values. Three 
methods have been investigated for this task: spectrum shifting, Higham remediation, and eigen- 
value clipping. When applied to the archived conjunctions involving NPD ECI state covariances, 
these three methods produced 2D P. values that differed only slightly from one another. The eigen- 
value clipping method has the advantages of being the simplest among the three, and the only one 
that provides a means of establishing reasonable, physics-based levels of remediation for NPD po- 
sition covariances. 
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